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TECHNICAL MEMORANDUM 


THE ANALYTICAL REPRESENTATION OF VISCOELASTIC MATERIAL 
PROPERTIES USING OPTIMIZATION TECHNIQUES 

INTRODUCTION 


One of the purposes of performing tests — either static, dynamic, mechanical, or electrical — is to 
determine the response of a system to a set of given test inputs. This test information can be reduced to a 
function of the system response, the dependent variable, to the test inputs, or the independent variables. 
This model can then be used to predict intermediate cases that are uneconomical to test. In the specific 
case of modeling the time-dependent properties of viscoelastic materials, it is particularly important that 
an accurate characterization of these properties be developed so that additional sources of error are not 
introduced into the solutions of problems involving these properties. 

A popular method used for formulating an analytical model for these time-dependent properties 
is to curve fit the test data to a predefined model, such as the Prony series. The Prony series can be 
defined as 


f(t) = A+IB i eK , (1) 

where t is defined as the independent variable, f(t) as the dependent variable, and A, B,-, and y are func- 
tion constants. The difficulty in reducing test data to the Prony series model is that the number of 
knowns in the model is exceeded by the number of unknowns. When this is encountered using tradi- 
tional methods, usually the only recourse is to assume values for some of the unknowns. 

In the case of the Prony series, the values for the exponential constants, y, are assumed to follow 
a decade pattern (i.e., y l =-l,y 2 = -10, y 3 = -100, ...). With the y terms defined, the constants A and B, 
can be solved for using least-squares regression techniques. The major drawback to this method is that 
the function is only as accurate as the values assumed for y. 

Rodriguez 1 presented a method in which the Prony series could be fit to test data using two 
Prony terms. This method allowed for the analytical solution of all constant terms. The only drawback to 
this method lay in the fact that the small number of Prony terms allowed for in this method may not 
characterize all the subtleties of the test data, thus the need for higher order terms in the function. 

The technique presented in this report allows data to be fit to the Prony series without the need 
for assuming values for the y terms. This technique determines the “optimum” Prony series function 
based on the input data and on minimizing the difference between the actual and calculated dependent 
variable. This report will expand on the method employed to determine the “optimum” function and will 
present some numerical examples to demonstrate the robustness, speed, and accuracy of the method. 

The computer program that was developed to perform the analysis was written in FORTRAN 77, 
and the analyses were performed on a 25-MHz 80386 PC with a math coprocessor installed. 



DESCRIPTION OF THE TECHNIQUE 


Recently, more and more emphasis has been placed on the design of structures that are light- 
weight yet offer high strength of stiffness. In order to achieve the minimum weight designs, optimization 
techniques and algorithms have been developed throughout the design community. One of the tools 
employed by the Structural Development Branch at Marshall Space Flight Center (MSFC) is the Design 
Optimization Tool (DOT) by VMA Engineering, Inc. 2 

DOT is a general optimization code written in FORTRAN 77 that will perform either constrained 
or unconstrained optimization, given the necessary design variables, constraints, and objective function. 
In order to access the optimization routines located within the DOT code, a header program is written 
that contains the definition of the initial DOT parameters; the initial value, lower bound, and upper 
bound for all design variables; a description of the objective function; a description of the functions 
required to constrain the design variables; and a calling statement to the DOT program. 

There are three optimization algorithms available with the DOT code for the solution of opti- 
mization problems. The first method, the Broydon-Fletcher-Goldfarb-Shanno (BFGS) algorithm, is 
employed only in unconstrained optimization problems. The other two available methods, the modified 
method of feasible directions (MMFD) and sequential linear programming (SLP), can be employed in 
the solution of constrained optimization problems. 

The technique developed in this report is an unconstrained optimization method and, therefore, 
employs the BFGS algorithm to minimize the objective function. The term “unconstrained optimization” 
means that the function being optimized (i.e., the objective function) is not restrained to any design 
space. Although this technique employs an unconstrained optimization method, it is not truly uncon- 
strained. The design variables (i.e., variables that can be manipulated in order to achieve a minimized or 
maximized objective function) are restrained through the use of side constraints which limit the range of 
values any one design variable can assume. The relation of constraints and side constraints to the objec- 
tive function are depicted in figure 1 (from ref. 2). 

In the evaluation of analytical data for an unrelated problem, it was found that the correlation of a 
function, derived through linear least-squares techniques, to a particular set of data could be improved 
through slight manipulation of the constant terms of the function. This led to the implementation of the 
DOT code to solve for the constants of a given function through the minimization of the sum of the 
residuals squared. The sum of the residuals squared gives a direct measurement of the accuracy of the 
curve fit and is defined as 

S r = lif(x t )-yi) 2 , ( 2 ) 

where /(x,) is defined as the value of the derived function at a given value of the independent variable, 
and y,- is the value of the dependent variable at the same independent variable value. BasicallyT the sum 
of the residuals squared reveals how close the function comes to passing through all the data points. This 
result is then used to calculate the correlation coefficient of the function. The correlation coefficient is 
also used to determine the accuracy of the curve fit, but reveals the percentage of error in the fit that is 
not explained by the function. The correlation coefficient is given by 



where S t is defined as the variance of the data about its mean. 


2 



x 2 



Figure 1. Relation of constraints to an objective function. 

Realizing the potential for this technique in curve fitting data to a variety of functions, it was 
decided to evaluate the technique by using it to reduce data to a function involving the Prony series. The 
Prony series provides a unique opportunity to evaluate the robustness of this technique because of the 
large number of design variables that can be applied against a given data set. For instance, to fit data to a 
Prony series with seven Prony terms, the optimization problem becomes an unconstrained minimization 
of a highly nonlinear function containing 15 design variables. Since the Prony series has two unknowns 
for each Prony term, the fi, and y terms, plus the constant term, A, the number of design variables can be 
determined by the relation 


Number of Design Variables = 2n+l , (4) 

where n is the number of Prony terms to fit to the data. 

Next, the issue of side constraints on the design variables will be discussed. In any regression of 
data, the constants of the function being determined have direct physical meaning to the independent and 
dependent variables. In the case of the Prony series being used to regress the time-dependent properties 
of viscoelastic materials, the A, and y, terms have a direct correspondence to the viscoelastic material 
properties. 

Evaluation of the physical significance of the constants requires examination of equation (1). In 
order for the exponential term in the series to be dimensionless, the y term will be required to have units 
of 1/time which is the form of a rate. Since the relaxation modulus of viscoelastic materials is a decaying 
function with respect to time, the y term is the rate of decay of the modulus, and thus will require a side 
constraint to ensure that it will always be negative. 

In equation (1), if time is allowed to approach infinity, then the series component of the function 
approaches zero and results in the constant A being the asymptotic value of the relaxation modulus at 
time infinity. Thus, in order to have physical significance, the constant A must always be positive. 
Therefore, a lower limit bound (i.e., a side constraint) of zero is set for the values that the constant A can 
assume. 
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Now, if time is allowed to approach zero, then all the exponential terms in the series component 
will approach one. Thus, the relaxation modulus at time zero is the summation of the constant A and the 
Bi terms. Therefore, the B, terms are also directly related to the relaxation modulus and, thus, require a 
side constraint on the lower bound to maintain positive values for all the B,- terms. 


COMPUTER IMPLEMENTATION 


Appendix A contains the FORTRAN 77 source code for the program PRONY which reduces 
data to the Prony series form based on the implementation of the DOT optimization subroutines. 

The program first initializes the variables necessary for the DOT subroutines to perform the 
optimization. The program then prompts the user for a data input filename, an output filename, and the 
number of Prony terms to be used. 

Next, the side constraints for each of the design variables are set to zero and either a large posi- 
tive or negative number, dependent upon whether the variable is allowed only positive or negative 
values, respectively. Then a call is made to a subroutine where the input file is opened, the xi and y t data 
pairs are read, and the mean value and variance of the dependent variable are calculated. Next, the loop 
is entered that controls the iterative process of performing the optimization task. Upon meeting the con- 
vergence criteria, the loop is exited, and the subroutine EVAL is evaluated to ensure that the objective 
function and the sum of the residuals squared are evaluated with the final values of the design variables. 
Finally, the output subroutine is called to send the results of the optimization to disk file. 

Now that the subject of design variables, constraints, side constraints, and the objective function 
have been discussed, it is time to briefly mention some of the other parameters that are required to per- 
form the optimization task through DOT. The first item to discuss is the convergence criteria. DOT uses 
two user-definable variables to control the convergence, the first being an absolute criteria and the 
second a relative criteria. 

If the initial value of the objective function is changed by less than the value for the absolute cri- 
teria during the optimization process, then the optimization task is assumed to have converged and will 
terminate. 2 The relative criteria will stop optimization when the relative change between the objective 
function in any two successive iterations is less than or equal to this variable. Since the objective func- 
tion generally assumes a large value on the initial calculation, because the initial values of the design 
variables are usually far from the optimum values, the Prony ser**' program relies on the relative criteria 
to control convergence. ^ 

The other remaining item to discuss concerning the implementation of the DOT code in the 
reduction of data to a Prony series function is the variable that controls the number of consecutive itera- 
tions that must satisfy the convergence criteria before optimization is terminated. By default, the DOT 
code sets this variable equal to two, since it is common to make little progress on one iteration, only to 
make significant progress on the next. It is suggested that if consistent progress is made on optimizing 
the objective function and if computer time is inexpensive, then this value can be raised. The PRONY 
program automatically overrides the default value and sets the variable equal to three, based upon test 
runs that showed additional improvement with the increased number of iterations. 
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NUMERICAL EXAMPLES 


Solid Rocket Motor TP-H1148 Propellant 

Modeling the relaxation modulus of solid rocket motor (SRM) TP-H1 148 propellant was 
selected as the first case since it was used as a numerical example by Rodriguez, 1 thus giving a direct 
comparison between the two methods. The data used by Rodriguez were obtained from a stress relax- 
ation test of SRM TP-H1 148 propellant that was tested at -30 °F and 2-percent strain. 3 The time varying 
relaxation moduli are listed in table 1. 

Table 1. Stress relaxation results of SRM TP-H1 148 propellant 
at -30 °F and 2-percent strain. 


Time (min) 

E r (lb/in 2 ) 

1.66X10- 3 

21,114 

6.61x10-3 

12,913 

1.66X10- 2 

10,992 

4.17x10-2 

10,114 

8.32x10-2 

8,399 

1.26x10-2 

7,566 

1.66X10- 1 

7,088 

1.66x10° 

3,568 

1.66X10 1 

1,981 


Rodriguez plotted the data points and fit a third-order B-spline curve to the data using an 
Intergraph computer-aided design (CAD) package called IGDS. Nineteen evenly spaced data points 
were then selected from the portion of the B-spline curve between 0.002 and 0.020 min. (Note: One of 
the requirements of the technique developed by Prony and used by Rodriguez is that the data points be 
spaced at equal time increments.) The constants A, B\, /i, Bi, and ji were determined to be 10,973.3, 
23,397.2, -380.8698, -79,056.7, and -1,971.8239, respectively. 

Comparing this curve-fitted Prony function to the actual test data through equations (2) and (3) 
yields a sum of residuals squared of 1.71xl0 8 and a correlation coefficient of 0.566609. On initial 
inspection, this method does not appear to generate a precise fit of the data points, but if the function is 
compared to the data points that basically fall within the aforementioned 0.002- and 0.020-min time 
interval, then a sum of residuals squared of 1.1 8x1 0 6 and a more impressive correlation coefficient of 
0.992128 are yielded. 

As a comparison of the Rodriguez method and the optimization method, a series of runs was 
made using the Prony program and the TP-H1 148 test data. These runs were made in an effort to deter- 
mine the sensitivity of the method to the initial starting point of the design variables and to the number 
of Prony terms utilized in the analysis. The output file listing generated by the Prony program for each 
of the cases listed in table 2 is contained within appendix C. Each file contains the date and time of the 
run, the input parameters, the sum of the residuals squared, the values for the A, B„ and ji terms, the 
correlation coefficient, and the central processing unit (CPU) time required for the analysis. 
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Table 2. Sensitivity study results for SRM TP-H1 148 propellant. 


Prony Terms 

Sums of Residuals 
Squared 

Correlation 

Coefficient 

Run Time (s) 

2 

1,427,300.1 

0.99716 

1.4 

3 

82,166.3 

0.99984 

3.0 

4 

82,201.7 

0.99984 

3.5 

5 

83,250.0 

0.99983 

6.6 


Referring to table 2, it can be seen that the first case of two Prony terms yielded a sum of residu- 
als squared that was 99.17 percent lower than that yielded by the Rodriguez method. The remaining 
three cases of three, four, and five Prony terms appear to have reached a lower bound on the sum of 
residuals squared that is approximately 94.24 percent lower than the case with two Prony terms. In all 
four cases analyzed, the CPU time required to run the analyses was insignificant due to the small number 
of data points available for the curve fit. 

Figures 2 and 3 graphically depict the relationship between the data points, the second-order 
Prony function derived from the Rodriguez method, and the second- and fourth-order Prony functions 
derived through the optimization technique, respectively. The figures show that the Rodriguez Prony 
function accurately models the propellant behavior over the time interval in which it was developed 
(0.002 to 0.020 min). It does not accurately model the late time behavior of the propellant. The second- 
and fourth-order optimized Prony functions, on the other hand, appear to accurately model the subtle 
fluctuations in the test data over the full-time interval. The apparent faceted nature of the curves in 
figures 2 and 3 is a direct result of the respective functions being plotted only at the test data points. 


V747-75 Viton Fluoroelastomer 

The second example also involves the time-dependent relaxation modulus of a viscoelastic 
material. This time the data were obtained from a stress relaxation test of the V747-75 Viton fluoro- 
elastomer. The stress relaxation test was conducted by the Polymers and Composites Branch at MSFC 
on an RDS-7700 dynamic spectrometer at 1-percent strain 4 and was shared with the author for the 
purpose of evaluating the subject technique. The data generated from this test and then used by the 
PRONY program are listed in appendix B. 

The data were originally fit to a 14-term Prony function that was developed in the classical 
method of assuming values for the ji terms, linearizing the resulting function, and then solving for the A 
and B t terms through least-squares techniques. The values for the constants A, and y terms *" 
determined from this technique are given in table 3. 

The resultant Prony function generated from applying these classically derived coefficients 
yields a sum of residuals squared of 9.20x1 0 6 and a correlation coefficient of 0.968443 when compared 
to the test data. This high correlation coefficient shows that the regressed function compares very well 
with the test data. 
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Table 3. Classically derived Prony coefficients. 


A: 

B i: 

111.800 

-45,830.000 

YT. -l.OxlO 7 

B 2 : 

2,547.000 

YT. -l.OxlO 6 

B 3 : 

3,156.000 

Yi'. -1.0x10 s 

B 4 : 

1,378.000 

y 4 : -l.OxlO 4 

B 5 : 

506.000 

75 : -l.OxlO 3 

B 6 : 

134.900 

/ 6 : -l.OxlO 2 

bt. 

86.450 

77 : -l.OxlO 1 

B 8 : 

47.140 

- 1 . 0 x 10 ° 

B 9 : 

38.420 

)§: -l.OxlO - 1 

# 10 : 

26.420 

7 i 0 : -l.OxlO - 2 

Bn'. 

24.340 

711 : -l.OxlO - 3 

B i 2 : 

11.540 

712 : -l.OxlO ^ 1 

B 13 : 

6.586 

7 i 3 : -l.OxlO - 5 

B 14 : 

15.010 

714 : -l.OxlO " 6 


A sensitivity study was again performed employing the PRONY program and the V747-75 Viton 
test data. The results of the study are contained in table 4. The output file listing generated by the 
PRONY program for each of the cases listed in table 4 is contained within appendix C. Each file 
contains the date and time of the run, the input parameters, the sum of the residuals squared, the values 
for the A, B t , and 7 terms, the correlation coefficient, and the CPU time required for the analysis. 


Table 4. Sensitivity study results for V747-75 Viton. 


Prony Terms 

Sums of Residuals 
Squared 

Correlation 

Coefficient 

Run Time (s) 

3 

478,677.4 

0.99838 

34.6 

4 

187,225.0 

0.99937 

64.4 

5 

58,402.1 

0.99980 

117.4 

6 

32,888.1 

0.99989 

205.0 

7 

19,110.1 

0.99994 

355.6 

8 

33,641.3 

0.99989 

334.3 

9 

139,630.9 

0.99953 

BHH19K 

10 

19,703.0 

0.99994 

661.7 

12 

19,705.8 

0.99993 

842.2 

14 

17,484.3 

0.99994 

1,330.7 


Referring to table 4, it can be seen that there is an indirect relationship between the number of 
terms and the sum of residuals squared up to and including the function containing seven Prony terms. 
Also, there exists a direct relationship between the number of terms and computer run time. With the 
exception of the nine Prony term case, it appears that after seven Prony terms, the optimized regression 
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has basically reached its lower limit for the sum of residuals squared. The small decreases in the sum of 
the residuals squared are not worth the additional run time of the program, but the CPU time on a PC is 
virtually free, and, if the user has the patience, the additional accuracy afforded with the 14th order 
Prony function may be worth the wait. An additional consideration that will affect the selection of the 
order of curve to use is the intended use of the function. If the function is to be used in “simple” hand 
calculations or short computer programs, then the higher-order curve may be selected. However, if the 
function is to be used in nonlinear stress analysis applications (i.e., applications in which a significant 
number of calculations will be performed with the function), one of the lower-order functions should be 
selected to reduce analysis times. 

Figures 4 and 5 graphically depict the relationship between the data points, the classically 
derived 14th order Prony function, and the 7th and 14th order Prony functions derived through the 
optimization technique, respectively. These two figures reveal the excellent correlation between the data 
and the optimally determined Prony functions. This correlation is most pronounced at the earlier time 
intervals since the higher relaxation moduli at these times artificially weight the optimization method 
and result in a better fit. At the later times, the relaxation moduli are smaller, and, thus, differences 
between the function and the data at these time intervals contribute less to the sum of the residuals 
squared (equation (2)). Therefore, at these later times, the function has more of a tendency to “wander” 
from the actual data than at the earlier times. 

CONCLUSIONS 

Application of optimization techniques to the task of reducing experimental data to an equivalent 
function of the form of the Prony series has been demonstrated. This technique generates functions that 
show excellent correlation to their original experimental data, as shown in the previous two numerical 
examples. It appears from the numerical examples that the technique is fairly robust, since it reduced 
data sets as small as 18 data pairs and as large as 99 data pairs. Also, the technique displayed the 
capability to reduce data sets containing data pairs spaced at varying time increments. 

It has been previously noted that this regression through optimization technique can generate 
functions that are artificially weighted to fit larger ordinate values more accurately than smaller ones. 
This artificial weighting was most apparent in the V747-75 Viton numerical example in which the 
function more closely approximated the experimental relaxation modulus at the early time intervals, 
where the modulus varied between 6,000 and 1,000 lb/in 2 , than at later times, where the modulus was 
less than 1,000 lb/in 2 . A few techniques were employed to attempt to minimize this residual weighting 
effect. Both normalization of the sum of residuals squared by the mean and changing the objective 
function so that the sum of the residuals would be minimized were attempted and generated poorer 
correlation of the resulting function to the experimental data. One method yet to be attempted is to 
decrease the frequency at which the large ordinate data pairs appear in the data set. Bower 5 states that 
this technique shows promise in a similar technique. ~ 

As was stated previously, the reduction of experimental data to accurate mathematical models is 
a vital part of most scientific and engineering research. This technique of regression through optimiza- 
tion can be applied to other mathematical models that are difficult to fit to experimental data through 
traditional regression techniques. The author has also employed this technique to fit an n independent 
variable power function to experimental data with a resultant function that more accurately predicts 
values for the dependent variable than a function generated through traditional regression techniques. In 
fact, the FORTRAN code for the program PRONY that is listed in appendix A would require minimal 
changes to determine coefficients for a wide range of mathematical models. 
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Figure 4. Comparison of Prony series functions, classical versus optimized methods. 



Figure 5. Comparison of Prony series functions, classical versus optimized methods. 
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APPENDIX A 


PRONY Source Code Listing 


C ****** Prony Curve Fit ****** 

C * Written By: Scott A. Hill * 

C ***************************** 

PARAMETER (converge-O. 00001 ) 

PARAMETER (itermax - 999) 

PARAMETER (itercon - 3) 

DOUBLE PRECISION r2,sumt,obj 

CHARACTER* 12 In, Out 

DIMENSION x(41),xl(41),xu(41) 

DIMENSION a (2000) ,y (2000) 

DIMENSION wk (2200) 

DIMENSION rprm (20) 

DIMENSION g (1) 

DIMENSION iwk(120) 

DIMENSION iprm (20) 

C INITIALIZE REAL AND INTEGER WORKSPACE FOR DOT 

DO 10 i-1,20 
rprm (i) =0.0 
iprm (i) -0 
10 CONTINUE 

rprm (3) - converge 
rprm (4) - converge 
iprm (3) - itermax 
iprm (4) - itercon 

C DEFINE # OF DESIGN VARIABLES, # OF CONSTRAINTS, AMT OF 
C DOT OUTPUT & WHETHER PROBLEM IS MINIMIZATION OR 
C MAXIMIZATION. 

nrwk-2200 

nriwk=120 

ncon=0 

iprint-0 

minmax— 1 

method-1 

C . DEFINE INITIAL VALUES AND UPPER/LOWER BOUNDS FOR DESIGN 
C VARIABLES 

WRITE (*,15) 

READ (*,'(A)') In 

15 FORMAT (' Input File: '\) 

WRITE (*,16) 

READ <*, ' (A)' ) Out 

16 FORMAT (' Output File: '\) 

WRITE (*,17) 

READ {*,*) ndv 

17 FORMAT ( # # of Prony Terms: # \) 

ndv « 2* ndv + 1 
x(l) - 150.0 
xl(l) - 0.0 
xu ( 1 ) - 1.0E+25 

DO 20 1-2, ndv, 2 
x (1 ) - 500.0 
x(i+l) - -3000.0 
20 CONTINUE 

DO 25 i - 2, ndv 
xl(i) - 0.0 
xu (i ) « 1.0E+25 

25 CONTINUE 


DO 30 i-3, ndv, 2 

xl(i) - -1.0E+25 
xu (i ) - 0.0 
30 CONTINUE 


C READ IN DATA & PERFORM STATISTICAL FUNCTIONS 
CALL InData (a,y, len, sumt) 

C READY TO OPTIMIZE 
info-0 

DO WHILE (.TRUE.) 


CALL DOT (info, method, iprint, iprint, ndv, 

. neon, x, xl, xu, ob j,minmax, g, rprm, 

. iprm,wk,nrwk,iwk,nriwk) 

CALL EVAL (ob j, x, y , a, sumt, len, r2, ndv) 

IF (info. EQ. 0) EXIT 

END DO 

C OPTIMIZATION COMPLETE — OUTPUT RESULTS 
CALL EVAL (ob j,x,y,a, sumt, len, r2, ndv) 

CALL OutPut (obj,x, ndv, r2, iprm (19) ,Type) 

END 

C* ************************************************** 

C SUBROUTINE TO EVALUATE FUNCTION AND CONSTRAINTS 

C* ************************************************** 

SUBROUTINE EVAL (ob j, x,y, a, sumt, len, r2, n, dv) 

REAL* 8 sumt, r 2, Value 

REAL* 8 ob j, fx 

DIMENSION x (* > 

DIMENSION y (*) ,a (*) 

obj - 0.0 

DO 100 i - 1, len 

fx - Value (x, a (i) , ndv) 

obj - obj + (y (i) - fx) * (y (i) - fx) 

100 CONTINUE 

r2 » (sumt - obj) /sumt 

RETURN 

END 

C* ******************************************** 

C FUNCTION TO CALCULATE VALUE OF PRONY SERIES 

C* ******** ******** **************************** 

REAL* 8 FUNCTION Value (array, x, n) 

REAL* 8 dv 

DIMENSION array (*) "" 

Value « array (1) 

DO 200 i - 2, n-1, 2 
dv - array (i+l)*x 
Value - Value + array (i) *DEXP (dv) 

200 CONTINUE 

END 
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Q* *************************************** 

C SUBROUTINE TO GET DATA FROM INPUT FILE 
£**************************************** 

SUBROUTINE InData {a, y, len, sumt > 

REAL* 8 mean, sumt 

CHARACTER* 12 In 

DIMENSION a (*) ,y (*) 

io - 0 
len - 0 
sum = 0.0 
sumt - 0.0 

OPEN (1 , FILE-In, STATUS-' OLD' ) 

DO 500 i - 1, 1000 

READ (1,*, END-510) a(i),y(±) 

500 CONTINUE 

510 CLOSE (1, STATUS-' KEEP' ) 

len - i - 1 

DO 520 i - l,len 
sum « sum + y(i) 

520 CONTINUE 

mean - sum/len 

DO 530 i - 1, len 

sumt - sumt + (y(i) - mean)*(y(i) - mean) 
530 CONTINUE 

RETURN 

END 


613 FORMAT (/5X, ' Prony Equation Form: y - A + ',A1, 

' Bxe A (' , Al, ' t ) ' ) 

620 FORMAT <//5X, ' Sura of the Residuals Squared: ', 
F14.4//',F14.4//12X,'A: ',F14.6> 

630 FORMAT <10X,'B', 12, ' : ' , F14 . 6, 15X, Al, 12, ' : ',F14.6) 
640 FORMAT </5X, 'Correlation Coefficient: ',F9.7// 

. 5X, ' # of Iterations 13) 

CLOSE (2, STATUS-' KEEP' ) 

RETURN 

END 


C* *************************************** 

C SUBROUTINE TO SEND DATA TO OUTPUT FILE 

C* *************************************** 

SUBROUTINE OutPut (ob j, x, ndv, r2, iter) 

REAL*8 obj,r2 

CHARACTER* 12 In, Out 

DIMENSION X (* ) 

OPEN (2, FILE-Out, STATUS-' UNKNOWN' ) 
k - 0 

WRITE (2,610) CHAR (228) , CHAR (231) 

WRITE (2,611) In 

WRITE (2,612) 150.0, 500. 0,CHAR(231), -1000.0, 
CHAR(231) ,0.0,0.000001 
WRITE (2, 620) ob j, x (1 ) 

DO 600 i - 2, ndv-1, 2 

k=k+l 

WRITE (2,630) k, x (i ), CHAR (231 ), k, x (i+1 ) 

600 CONTINUE 

WRITE (2, 640) DSQRT (r2 ) , iter 

610 FORMAT (/5X, ' Prony Equation Form: y - A + ', 

Al, ' Be A (' , Al , ' t) ' ) 

611 FORMAT (/5X,' Input File: ',A) 

612 FORMAT (//5X, ' Initial Value for A: ',F12.4/5X, 

'Initial Value for B: ', F12 . 4/5X, ' Initial Value', 
. ' for ' , A, ' : ', F12 . 4/5X, ' Upper Limit for ',A, 

. ': ' , F12 . 4 //5X, ' Convergence Criteria: ' 

, F12 .10) 
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APPENDIX B 


V747-75 Vlton Stress Relaxation Data 


Time (min) 

E r- (PSj) 

Time (min) 

ir (psi). 

Time (min) 

3 . 6700E-06 

5 . 6500E+03 

2 . 0650E+00 

2 . 4780E+02 

1.1620E+06 

4 . 8090E-06 

5 . 3240E+03 

2 . 7050E+00 

2 . 4250E+02 

1 . 5220E+06 

6.3000E-06 

4 . 8970E+03 

3 . 5450E+00 

2 . 3740E+02 


8.2550E-06 

4 . 4720E+03 

4 . 6450E+00 

2 . 3130E+02 


1 . 0820E-05 

4.1130E+03 

6 . 0850E+00 

2 . 2750E+02 


1 . 4170E-05 

3 . 7880E+03 

7.9730E+00 

2 . 2220E+02 


1 . 8570E-05 

3 . 4710E+03 

1 . 0450E+01 

2. 1800E+02 


2 . 4330E-05 

3 . 1760E+03 

1 . 3690E+01 

2 . 1390E+02 


3.1880E-05 

2 . 8910E+03 

1.7930E+01 

2 . 0900E+02 


4.1770E-05 

2 . 6240E+03 

2 . 3500E+01 

2 . 0430E+02 


5 . 4720E-05 

2 . 3760E+03 

3 . 0790E+01 

2 . 0000E+02 


7 . 1700E-05 

2 . 1530E+03 

4 . 0340E+01 

1 . 9600E+02 


9 . 3950E-05 

1 . 9470E+03 

5 . 2860E+01 

1 . 9210E+02 


1 . 2310E-04 

1 . 7710E+03 

6. 9260E+01 

1 . 8770E+02 


1 . 6130E-04 

1 . 6130E+03 

9 . 0740E+01 

1 . 8630E+02 


2 . 1130E-04 

1 . 4810E+03 

1 . 1890E+02 

1 . 8490E+02 


2 . 7 690E-04 

1 . 3530E+03 

1.5580E+02 

1 . 8310E+02 


3 . 6280E-04 

1 . 1850E+03 

2 . 0410E+02 

1.7850E+02 


4 . 7530E-04 

1 . 0590E+03 

2.6740E+02 

1 . 7390E+02 


6 . 2280E-04 

9 . 4020E+02 

3 . 5040E+02 

1 . 7080E+02 


8 . 1600E-04 

8 . 5680E+02 

4 . 5910E+02 

1 . 6780E+02 


1 . 0690E-03 

7 . 8990E+02 

6 . 0150E+02 

1 . 6500E+02 


1.4010E-03 

7 . 3480E+02 

7 . 8820E+02 

1 . 6060E+02 


1 . 8350E-03 

6 . 8780E+02 

1 . 0330E+03 

1 . 5820E+02 


2 . 4050E-03 

6 . 4450E+02 

1 . 3530E+03 

1 . 5670E+02 


3 . 1510E-03 

6 . 0100E+02 

1 . 7730E+03 

1 . 5310E+02 


4 . 1290E-03 

5 . 5750E+02 

2 . 3230E+03 

1 . 5140E+02 


5.4090E-03 

5 . 2570E+02 

3 . 0440E+03 

1 . 4950E+02 


7 . 0880E-03 

4 . 9600E+02 

3 . 9880E+03 

1 . 4680E+02 


9 . 2870E-03 

4 . 6860E+02 

5.2250E+03 

1 . 4510E+02 


1 . 2170E-02 

4 . 4330E+02 

6 . 8460E+03 

1 . 4330E+02 


1 . 5940E-02 

4 . 2210E+02 

8 . 9700E+03 

1 . 4160E+02 


2 . 0890E-02 

4 . 0390E+02 

1 . 1750E+04 

1 . 3980E+02 


2 . 7370E-02 

3 . 8690E+02 

1 . 5400E+04 

1 . 3810E+02 


3 . 5860E-02 

3 . 7290E+02 

2 . 0180E+04 

1 . 3650E+02 


4 . 6990E-02 

3. 6500E+02 

2 . 6440E+04 

1 . 3600E+02 


6 . 1560E-02 

3.5740E+02 

3 . 4640E+04 

1 . 3450E+02 


8 . 0660E-02 

3 . 4750E+02 

4 . 5380E+04 

1.3350E+02 


1 . 0570E-01 

3 . 3350E+02 

5. 9460E+04 

1 . 3230E+02 


1 . 3850E-01 

3 . 2100E+02 

7 . 7910E+04 

1 . 3130E+02 


1.8140E-01 

3 . 1100E+02 

1 . 0210E+05 

1 . 2990E+02 


2 . 3770E-01 

3 . 0070E+02 

1 . 3380E+05 

1 . 2910E+02 


3 . 1150E-01 

2 . 9330E+02 

1 . 7520E+05 

1 . 2830E+02 


4 . 0810E-01 

2.8600E+02 

2 . 2960E+05 

1.2770E+02 


5.3470E-01 

2 . 7680E+02 

3 . 0090E+05 

1 .2660E+02 


7 . 0060E-01 

2 . 7220E+02 

3 . 9420E+05 

1 . 2540E+02 


9 . 1800E-01 

2 . 6610E+02 

5 . 1650E+05 

1 . 2430E+02 


1 .2030E+00 

2 . 5990E+02 

6 . 7 670E+05 

1 . 2320E+02 


1 . 5760E+00 

2 . 5340E+02 

8 . 8670E+05 

1.2160E+02 



E - (P3i) 

1 .2050E+02 
1 . 1900E+02 
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APPENDIX C 


PRONY Program Output Files 
TP-H1148 SRM Propellant 


07/28/1992 07/28/1992 

10:35:48 10:49:23 


Prony Equation Form: y 

Input File: srm2.dat 


T xt \ 

r 

{- 1 

Prony Equation Form: y - A + £ Be 1 


Input File: srm2.dat 


Initial Value for A: 
Initial Value for B: 
Initial Value for T: 
Upper Limit for t: 


2700.0000 

8000.0000 
-300.0000 

.0000 


Initial Value for A: 
Initial Value for B: 
Initial Value for t: 
Upper Limit for t: 


100.0000 

100.0000 

- 100.0000 

.0000 


Convergence Criteria: .0000000100 


Convergence Criteria: .0000010000 


Sum of the Residuals Squared: 

1427300.1201 

Sum of 

the Residuals Squared: 

82201.7363 

A: 2779.583000 



A: 

1982.994000 



Bl: 8 . 794046E+03 

Tl: 

“4.731266E+00 

Bl: 

5.718736E+03 

Tl : 

-9. 948862E+00 

B2: 1 . 74802 9E+04 

x2: 

-3.604886E+02 

B2 : 

4 . 273885E+03 

t2 : 

-5.971107E-01 




B3: 

7 . 62974 5E+03 

t3: 

-3 . 896799E+02 

Correlation Coefficient: 

.9971610 

B4 : 

1 . 012811E+04 

t4 : 

-3 . 969893E+02 

# of Iterations: 36 



Correlation Coefficient: 

.9998367 

Optimization Time: Ohr 

Omin 

1.4sec 








# of Iterations: 45 






Optimization Time: Ohr 

Omin 

3. 5sec 


07/28/1992 


07/28/1992 


10:44:01 


11:09:14 



r 



Prony Equation Form: 

y «* A + £ Be 1 

Prony Equation Form: 

y « A + £^Be 

Input File: srm2.dat 


Input File: srm2.dat 


Initial Value for A: 

1000.0000 

Initial Value for A: 

1400.0000 

Initial Value for B: 

100.0000 

Initial Value for B: 

131.0000 

Initial Value for t: 

“100.0000 

Initial Value for T: 

-100.0000 

Upper Limit for t: 

.0000 

Upper Limit for t: 

.0000 

Convergence Criteria: 

.0000010000 

Convergence Criteria: 

.0000010000 


Sum of the Residuals Squared: 82166.2936 Sum of the Residuals Squared: 83249.9692 


A: 

1980.808000 



Bl: 

1.775723E+04 

xl: 

-3. 946309E+02 

B2: 

4 . 318967E+03 

t2: 

-6. 028311E-01 

B3: 

5 . 687424E+03 

t3: 

-1 . 008739E+01 

Correlation Coefficient: 

.9998368 


A: 

1984.816000 



Bl: 

3. 686836E+03 

Tl: 

-I . 063451E+01 

B2: 

2.206965E+03 

x2 : 

-7.734952E+00 

B3: 

4 . 170810E-01 

t3 : 

-2 . 230149E+03 

B4: 

4 . 079599E+03 

t4: 

-5 . 694004E-01 

B5: 

1 . 772458E+04 

t5 : 

-3.917794E+02 


# of Iterations: 54 Correlation Coefficient: .9998346 

Optimization Time: Ohr Omin 3.0sec 

# of Iterations: 65 

Optimization Time: Ohr Omin 6.6sec 
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V747-75 Viton Fluoroelastomer 


01/22/1992 

16:12:15 


Prony Equation Form: y 


A + 



Input File: data.rlx 


01/23/1992 

07:49:44 


Prony Equation Form: y 


A + 



Input File: data.rlx 


Initial Value for A: 150.0000 
Initial Value for Bx 500.0000 
Initial Value for x: -500.0000 
Upper Limit for x: .0000 

Convergence Criteria: .0000100000 


Initial Value for A: 150.0000 
Initial Value for B: 500.0000 
Initial Value for x: -1000.0000 
Upper Limit for x: .0000 

Convergence Criteria: .0001000000 


Sum of the Residuals Squared: 478677.4458 

A: 187.385400 

Bl: 3.819927E+03 xl: -7 . 623515E+04 

B2: 4 . 861179E+02 x2: -2.872654E+01 

B3: 1 . 974550E+03 x3: -3. 894463E+03 

Correlation Coefficient: .9983830 

# of Iterations: 93 

Optimization Time: Ohr Omin 34.6sec 


Sum of the Residuals Squared: 58402.0883 


A: 

150.559100 



Bl: 

3. 135792E+03 

xl: 

-1 . 681769E+05 

B2: 

3. 978183E+02 

x2 : 

-8 . 580777E+01 

B3: 

1 . 601320E+02 

X3: 

-1.008694E-01 

B4: 

1 . 298270E+03 

x4 : 

-2 . 538760E+03 

B5: 

2. 147726E+03 

x5: 

-2 . 4 67912E+04 


Correlation Coefficient: .9998029 

♦ of Iterations: 154 

Optimization Time: Ohr lmin 57.4sec 


01/22/1992 

16:15:53 


Prony Equation Form: y 


A + 



Input File: data.rlx 


01/23/1992 

07:53:54 


Prony Equation Form: y 


A + 



Input File: data.rlx 


Initial Value for A: 150.0000 
Initial Value for B: 500.0000 
Initial Value for X: -500.0000 
Upper Limit for x: .0000 

Convergence Criteria: .0000100000 


Initial Value for A: 150.0000 
Initial Value for B: 500.0000 
Initial Value for x: -1000.0000 
Upper Limit for x: .0000 

Convergence Criteria: .0000100000 


Sum of the Residuals Squared: 187224.9701 

A: 159.108800 

Bl: 2 . 115458E+02 xl: -5. 073118E-01 

B2: 6 . 906992E+02 x2: -3.441239E+02 

B3: 2 . 04 6774E+03 x3: -7 . 607188E+03 

B4: 3 . 635521E+03 x4: -9.960206E+04 

Correlation Coefficient: .9993678 

# of Iterations: 114 

Optimization Time: Ohr lmin 4.4sec 


Sum of the Residuals Squared: 32888.0913 


A: 

136.880500 



Bl: 

8. 903401E+01 

Xl: 

-3. 697306E-03 

B2: 

3.241039E+03 

x2 : 

-1 . 343083E+05 

B3: 

1 . 139755E+03 

x3 : 

-2.481861E+03 

B4: 

1 . 876921E+03 

x4 : 

-1 . 908836E+04 

B5: 

3. 799151E+02 

x5: 

-1 . 5124 57E+02 

B6: 

1 . 4 91300E+02 

x6 : 

-2.062958E+00 


Correlation Coefficient: .9998890 

♦ of Iterations: 202 

Optimization Time: Ohr 3min 25.0sec 
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01/23/1992 

08:00:47 


Prony Equation Form: y 


A + £ 


K'l 


Input File: data.rlx 


Initial Value for A: 
Initial Value for B: 
Initial Value for x: 
Upper Limit for x: 


150.0000 

500.0000 
- 1000.0000 

.0000 


Convergence Criteria: .0000100000 


01/23/1992 

08:21:34 


Prony Equation Form: y 


A + I 


r Tt i 

r j 


Input File: data.rlx 


Initial Value for A: 
Initial Value for B: 
Initial Value for x: 
Upper Limit for x: 


150.0000 

500.0000 
- 1000.0000 

.0000 


Convergence Criteria: .0000100000 


Sum of the Residuals Squared: 19110.1326 


A: 

137.237000 



Bl: 

9.063905E+01 

Xl: 

-4 . 059385E-03 

B2 : 

1 . 2 96831E+03 

x2 : 

-2 . 984 362E+03 

B3: 

1 . 1 64508E+03 

x3: 

-2.645348E+04 

B4 : 

9 . 470085E+02 

x4 : 

-2 . 577721E+04 

B5: 

1 . 540934E+02 

x5: 

-2 . 360951E+00 

B6 : 

4 . 0704 35E+02 

x6: 

-1.741451E+02 

B7: 

3.099127E+03 

x7 : 

-1.704219E+05 


Correlation Coefficient: .9999355 

# of Iterations: 263 

Optimization Time: Ohr 5min 55.6sec 


Sum of the Residuals Squared: 139630.8670 


A: 

139.145000 



Bl: 

8. 233751E+01 

xl: 

-7 . 876419E-03 

B2: 

4 . 220577E+02 

t2: 

-1.708051E+02 

B3: 

1 . 0987 96E+03 

t3: 

-3 . 100827E+03 

B4: 

2 . 08604 9E+02 

t4 : 

-1 . 831882E+04 

B5: 

3 . 555962E+03 

x5: 

-7 . 973099E+04 

B6: 

3 . Q27768E+02 

x6: 

-1 . 972911E+04 

B7: 

1 . 4 17834E+02 

x7 ; 

-2 . 247524E+00 

B8: 

5 . Q97406E+02 

x8 : 

-5. 760359E+03 

B9: 

1 . 100006E+01 

x9: 

-5. 820372E-04 


Correlation Coefficient: .9995286 

# of Iterations: 215 

Optimization Time: Ohr 6min 50.6sec 


01/23/1992 

08:07:17 

Prony Equation Form: y 

Input File: data.rlx 


A + 



01/23/1992 

08:33:42 

Prony Equation Form: y 

Input File: data.rlx 


A + £ 



Initial Value for A: 
Initial Value for B: 
Initial Value for x: 
Upper Limit for x: 


150.0000 

500.0000 
- 1000.0000 

.0000 


Convergence Criteria: .0000100000 


Initial Value for A: 
Initial Value for B: 
Initial Value for x: 
Upper Limit for X: 


150.0000 

500.0000 
- 1000.0000 

.0000 


Convergence Criteria: .0000100000 


Sum of the Residuals Squared: 


19073.0373 


of 

the Residuals Squared: 

33641.3004 









A: 

132.907500 



A: 

137.107400 



Bl: 

6 . 29991 5E+01 

xl : 

-9 . 576948E-04 

Bl: 

1.521836E+02 

xl: 

-2.256213E+00 

B2 : 

8.757581E+01 

x2 : 

-1 . 771211E-01 

B2: 

1 . 4 99064E+03 

x2 : 

-2 . 030280E+04 

B3: 

5.747751E+02 

x3: 

-2 . 629937E+03 

B3: 

9 . 00594 7E+01 

x3 : 

-3. 925191E-03 

B4: 

1 . 4 20293E+02 

x4 : 

-9 . 763288E+00 

B4: 

3 . 887046E+02 

x4 : 

-1 . 622622E+02 

B5: 

6.023077E+02 

x5: 

-1.751084E+04 

B5: 

3 . 1 634 55E + 01 

x5: 

-1 . 420Q35E+04 

B6: 

3.735197E+02 

x6: 

-2 . 204888E+02 

B6: 

1 . 140985E+03 

x6: 

-2 . 5662 99E+03 

B7: 

5. 726660E+02 

x7 : 

-2.89067.6E+03 

B7: 

3 . 2 981 98E+02 

x7 : 

-1.568566E+04 

B8: 

1 . 24 3688E+03 

x8 : 

-2 . 652360E+04 

B8: 

3.214 954 E + 03 

x8 : 

-1 . 324 574 E+05 

B9: 

3.099371E+03 

x9 : 

-1 . 429671E+05 





B10: 

1 . 601 951E+02 

xlO : 

: -9. 859158E+03 

relation Coefficient: 

.9998864 









Correlation Coefficient: 

.9999356 

: Iterations: 207 







.mization Time: Ohr 

5min 

34 . 3sec 

# of Iterations: 289 







Optimization Time: Ohr 

llmin 

1 . 7sec 
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01/23/1992 

11:34:44 


01/23/1992 

12:23:32 


Prony Equation Form: y 


A + 



Input File: data.rlx 


Initial Value for A: 
Initial Value for B: 
Initial Value for x: 
Upper Limit for x: 


150.0000 

500.0000 
- 1000.0000 

.0000 


Convergence Criteria: .0000100000 


Sum of the Residuals Squared: 19705.8298 


A: 

Bl: 

133.527800 
1 . 482939E+02 

xl: 

-1.012032E+01 

B2: 

8 . 787807E+01 

x2 : 

—1 . 609448E-01 

B3: 

6 . 217831E+01 

x3: 

-8. 913591E-04 

B4: 

3. 951779E+02 

x4: 

— 2.483659E+02 

B5: 

2.622941E+02 

x5: 

-6.500080E+03 

B6: 

1 . 728592E+03 

x6: 

-2. 389782E+04 

B7: 

4.072714E+01 

x7 : 

-3. 761625E+03 

B8 : 

9 . 704579E+01 

x8: 

-9. 570082E+03 

B9: 

9 . 301317E+02 

x9: 

-2. 696142E+03 

B10: 

3. 821 681E+01 

xlO: 

-3. 602862E+03 

Bll: 

1 . 437028E-03 

xll: 

-2.693170E+04 

B12: 

3 . 131319E+03 

xl2 : 

-1 . 417388E+05 


Correlation Coefficient: .9999335 


Prony Equation Form: y 


A + 



Input File: data.rlx 


Initial Value for A: 
Initial Value for B: 
Initial Value for x: 
Upper Limit for x: 


150.0000 

500.0000 
-3000.0000 

.0000 


Convergence Criteria: .0000100000 


Sum of the Residuals Squared: 17484.2568 


A: 

136.978600 

Bl: 

7.829514E+01 

B2: 

1 . 521367E+02 

B3: 

3.104636E+03 

B4: 

1 . 371974E-01 

B5: 

6. 143808E+01 

B6; 

3. 986375E+02 

B7: 

2 . 540747E+01 

B8: 

1 . 882670E-03 

B9: 

8. 968185E+01 

B10: 

8. 744885E+00 

Bll: 

1.186864E+03 

B12: 

9. 729550E+02 

B13: 

1 . 141888E+03 

B14: 

7 . 512561E+01 


Correlation Coefficient: 


Xl: -6. 173314E+03 
X2 : -2 . 203388E+00 
x3: -1 . 868573E+05 
x4 : -1 . 328096E+04 
x5 : -7 . 706898E+03 
X6: -1.654526E+02 
x7: -4 . 293449E+03 
x8: -5.806276E+03 
x9 : -3. 824855E-03 
xlO: -1 . 269737E+04 
xll: -2 . 782989E+03 
xl2 : -2 . 831181E+04 
X13: — 2 . 85991 6E+04 
Xl4 : -1 . 803128E+04 

.9999410 


# of Iterations: 263 

Optimization Time: Ohr 14min 2.2sec 


# of Iterations: 322 

Optimization Time: Ohr 22min 10.7sec 
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